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The many-body formalism for dynamical mean-field theory is extended to treat nonequilibrium 
problems. We illustrate how the formalism works by examining the transient decay of the oscillating 
current that is driven by a large electric field turned on at time t — 0. We show how the Bloch 
oscillations are quenched by the electron-electron interactions, and how their character changes 
dramatically for a Mott insulator. 
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Introduction. Dynamical mean-field theory (DMFT) 
was introduced in 1989 as a technique to solve the quan- 
tum many-body problem by taking the limit where the 
number of spatial dimensions goes to infinity jlj ■ In this 
limit, with a proper scaling of the hopping matrix el- 
ements, the electron-electron correlations are described 
by a local self-energy. Hence the many-body problem on 
a lattice is mapped onto an effective many-body problem 
for a single-site impurity (in a time-dependent field), with 
a self-consistency condition that fixes the time-dependent 
field so that the Green's function for the impurity is 
identical with the local Green's function for the lattice. 
Since then, DMFT has been employed to solve virtually 
all many-body problems described by model Hamiltoni- 
ans Q, has been generalized to describe strong electron 
correlations in real materials |3( and to describe inhomo- 
geneous systems 0, Q . All of this work has focused on 
the equilibrium case. In this contribution, we illustrate 
how to generalize DMFT to nonequilibrium situations, 
and we present results for how the Bloch oscillations of 
a strongly correlated material are quenched by electron- 
electron interactions, and how their character changes 
after the Mott metal-insulator transition. 

Bloch 6] and Zener Q theorized that electrons on 
a lattice undergo an oscillatory motion when placed 
in a uniform static electric field, because the electron 
wavevector, which evolves under the electric field, is 
Bragg reflected whenever it reaches a Brillouin zone 
boundary. But in metals, Bloch oscillations have never 
been seen, because the electron relaxation time is so 
short, the electrons are scattered before they reach the 
zone boundary and Bragg reflect. Bloch oscillations have 
been observed in semiconducting heterostructures ||, 
Josephson junctions and cold-atom systems 

Formalism. The many-body formalism for nonequi- 
librium dynamical mean-field theory is straightfor- 
ward to dev elop within the Kadanoff-Baym-Keldysh ap- 
proach [ill Il2j |. Because nonequilibrium problems are 
not time-translation invariant, we need to work with 
Green's functions that depend on two times. There 
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FIG. 1: Kadanoff-Baym-Keldysh contour for the two-time 
Green's functions in the nonequilibrium case. We take the 
contour to run from —5 to t ma x and back, and then extends 
downward parallel to the imaginary axis for a distance of f3. 
The field is turned on at t = 0; i.e., the vector potential is 
nonzero only for positive times. 



are two independent Green's functions that need to be 
determined — the retarded Green's function G R , which 
describes the density of available quantum-mechanical 
states, and the lesser Green's function G < , which deter- 
mines how electrons occupy those quantum states. Both 
Green's functions can be extracted from the contour- 
ordered Green's function, which is defined for any two 
time values that lie on the Kadanoff-Baym-Keldysh con- 
tour shown in Fig. ^ We imagine our system to be 
in equilibrium until time t = where a field is turned 
on. The contour starts at some time before the field is 
turned on, runs out to a maximal time, then returns to 
the original time, and finally moves parallel to the nega- 
tive imaginary axis a distance (equal to the inverse of 
the temperature of the original equilibrium distribution). 

Since the many-body perturbation theory diagrams 
are identical in structure for equilibrium and nonequi- 
librium perturbation theories [13j , the perturbative anal- 
ysis of Metzner 0] guarantees that the nonequilibrium 
self-energy is also local in DMFT. Hence, the nonequilib- 
rium DMFT problem can be mapped onto an impurity 
problem in time-dependent fields, just like the equilib- 
rium problem, except now the fields have two time argu- 
ments. The basic structure of the iterative approach to 
solving the DMFT equations 0] continues to hold. We 
start with a guess for the self-energy (which is usually 
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chosen to be equal to the equilibrium self-energy), then 
we sum the momentum-dependent Green's function over 
the Brillouin zone to produce the local Green's function. 
Next the dynamical mean-field for the impurity prob- 
lem is extracted by using Dyson's equation for the local 
Green's function and self-energy, the impurity problem is 
solved in the dynamical mean-field to produce the impu- 
rity Green's function, and Dyson's equation is used again 
to extract the impurity self-energy. In the self-consistent 
solution of the DMFT equations, the impurity self-energy 
will be equal to the lattice self-energy. If they are differ- 
ent, then the new lattice self-energy is taken to be equal 
to the new impurity self-energy, and the loop is iterated 
until it converges. The nonequilibrium algorithm is mod- 
ified as follows: (i) the summation over the Brillouin zone 
now requires at least a double integral over two band en- 
ergies; (ii) the Green's functions are described by discrete 
matrices with time indices that run over the contour; and 
(iii) the impurity problem solver must be generalized to 
the nonequilibrium case. 

For concreteness, we assume the electric field E(t) is 
spatially uniform, but can depend on time (we assume 
that the magnetic field, however, is small, and neglect 
all magnetic-field effects). We work in the Hamiltonian 
gauge, where the scalar potential vanishes, and the elec- 
tric field is determined by a time derivative of the vec- 
tor potential E(£) = —dA(t)/dt, in units where c = 1. 
The noninteracting problem of Bloch electrons in an elec- 
tric field can be solved exactly by using the Peierls sub- 
stitution 0, 0], and if we take the electric field to 
lie along the diagonal direction, then the noninteracting 
momentum-dependent Green's functions on the lattice 
depend only on two explicit functions of momentum 

e k = -^gcosk i; £ - k = -^!>k 4 , (1) 

rather than all components of the momentum. Here we 
set the lattice constant a equal to 1, and we consider the 
case of nearest-neighbor hopping on a hypercubic lattice 
in d-dimensions with a hopping parameter t = t* /2y/d; t* 
will be taken as the energy unit. In the limit d — > oo, the 
two "band energies" are distributed with a joint Gaussian 
density of states 

p(ek,e k ) = - e - e k e -efc. (2) 

IT 

In the interacting case, the dressed contour-ordered 
Green's function satisfies Dyson's equation, with a local 
self-energy, so that 

G(k, t, t') = [Go 1 (k, t , t') ~ E(t, f )] _1 , (3) 

where the Green's functions and self-energy are continu- 
ous matrix operators defined on the contour (i. e., the 
time indices of the matrices run along the contour), and 



the —1 superscripts denote the matrix inverse of the re- 
spective operators. The noninteracting Green's function 
in a field (for both times larger than 0; one can easily 
work out the generalizations for cases when the field has 
not been turned on) is 

G (k, t, f) = i[/(e k - ,i) - 0&, OJe-W*-*'' (4) 

w — ie^fsin eEt — sin eEt')/eE 
X c 

— iekfcos eEt— cos eEt')/eE 
X 6 , 

where e is the electron charge, E is one component of 
the electric field along a Cartesian axis (all components 
are equal for a field directed along the diagonal), and 
f(x) = l/[l + cxp(x)] is the Fcrmi-Dirac distribution (we 
set h = 1). The symbol 8 c (t,t') is equal to one if t is 
ahead of t' on the contour and is zero otherwise. Cal- 
culating the local Green's function requires evaluating a 
two-dimensional integral over p(ek, e k ) of a matrix- valued 
integrand, which requires a matrix inversion to determine 
it. Once the local Green's function G has been found, 
we use Dyson's equation to extract the dynamical mean- 
field, denoted X(t,t r ), which satisfies 

A(t, = {idf + ii)5 c (t, - G-\t, t 1 ) - E(t, 0, (5) 

where the derivative with respect to time is taken along 
the contour, and the delta function is defined on the con- 
tour such that j c dt5 c (t, t')F(t) = F(t'). 

Next, the impurity problem must be solved for elec- 
trons evolving in the dynamical mean field. In general, al- 
gorithms have not yet been developed to solve this prob- 
lem for all Hamiltonians, but the impurity problem can 
be immediately solved for the spinless Falicov-Kimball 
model 19], which involves single-band conduction elec- 
trons hopping on a lattice, and localized electrons which 
do not move but do interact with the conduction elec- 
trons when they are in the same unit cell via a screened 
Coulomb interaction U. The Hamiltonian (in the absence 
of a field) is then 

n = -^rr d E( c h- + c h) + uJ2 w * c W (6) 

fa'} » 

Here, we have c\ (q) create (annihilate) a spinless con- 
duction electron at site i and w% — or 1 is the localized 
electron number operator at site i. Although the Falicov- 
Kimball model is a simple many-body physics model, 
it does have a Mott metal-insulator transition (but the 
model does not include Zcncr tunneling because there are 
no higher energy bands). The solution to the impurity 
problem can be found by solving the equations of motion 
for the contour-ordered Green's function resulting in 

G imp (t,t') = (i-t«i)[(i^+A*)*c(*,t')-A(*.or 1 

+ Wl [{id c t U)S c (t, t') - \{t, t')}- 1 (7) 
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with u>i the average localized electron filling. The Dyson 
equation in Eq. J5| is then employed to extract the im- 
purity self-energy, and the algorithm is iterated until it 
converges. 

There are a number of important technical details in 
these calculations. First, we discretize the contour (we 
choose a real-time spacing At which varies from 0.1 to 
0.0333, and we fix the spacing along the imaginary axis 
to At = O.li) and evaluate integrals over the contour 
by discrete summations using the midpoint rectangular 
integration rule. The matrix operators are general com- 
plex matrices, which are manipulated using standard lin- 
ear algebra packages (the largest matrices used are about 
2200 x 2200). In addition, the delta function changes sign 
along the negative real time branch, and is imaginary 
along the last branch of the contour while the deriva- 
tive of the delta function is evaluated by a two-point dis- 
cretization involving the diagonal and the first subdiago- 
nal, but one also needs to include one matrix element at 
the upper right corner to preserve the proper boundary 
conditions. 

We perform the two-dimensional energy integration by 
Gaussian quadrature with both 100 and 101 points in 
each dimension, and we average the two results. Since 
the calculation of each matrix in the integrand of the inte- 
gral is independent of every other quadrature point, this 
part of the code is easily parallelized. We require 20,201 
matrix inversions for each DMFT iteration. The impu- 
rity solver is a serial code, that cannot be parallelized, 
because the matrix operations need to all be performed 
in turn. We typically require between ten and fifty itera- 
tions to reach convergence of the results (the total com- 
puter time for the calculations presented here was about 
600,000 cpu-hours on a Cray XT3). Once converged, we 
calculate the current by evaluating the operator average 

=- e ^v(k- e M)G < (k,t,i). (8) 

k 

The velocity component is Vi (k) = t*sin(ki)/2Vd, and 
all components of the current are equal when the field lies 
along the diagonal. We also calculate the equal time re- 
tarded and lesser Green's functions and their first two 
derivatives and compare the results to the exact val- 
ues |2jj. In general, these "moments" are quite accurate 
as the step size is made smaller, and we find that if we 
use a Lagrange interpolation formula to extrapolate the 
results to At = 0, we can achieve even higher accuracy 
for most values of U. Details of these numerical issues 
and of the accuracies will be presented elsewhere 12 lj . 

Numerical results. We produce numerical calculations 
of the nonequilibrium current as a function of time for 
the case of half-filling, where the conduction electron and 
the localized electron fillings are each equal to 0.5. This 
system has a metal-insulator transition at U — \pl. In 
the case where there is no scattering (U — 0), the Bloch 
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FIG. 2: Scaled nonequilibrium current for different values 
of U. Note how the Bloch oscillations are rapidly reduced in 
amplitude as the scattering increases. (Color on-line.) 

oscillations continue forever. In the presence of scatter- 
ing, the Bloch oscillations maintain the same approxi- 
mate "periodicity", but the amplitude decays. In Fig. El 
we plot the current for the noninteracting case, the case 
of a strongly scattering metal (U — 0.5, red), the case 
of an anomalous metal (U = 1, green), and of a near 
critical insulator (U — 1.5, blue). The initial temper- 
ature of the system satisfied j3 — 10, and the field is 
turned on at t — 0. The electric field is set equal to 
one in magnitude, E = 1. Note how the Bloch oscilla- 
tions are damped as the scattering increases. Although 
a Boltmann equation approach always predicts that the 
oscillations are damped and disappear on a time-scale on 
the order of the relaxation time, and approach a constant 
steady state, we do not see this full evolution within the 
time-window that we performed these calculations. Most 
of the data given here involve a scaling of the data with 
At = 0.1, 0.067, 0.05, and 0.04 to the At -> limit. 
Note that the data for a fixed step size in time always 
shows a small current for t < 0, but when scaled, the 
current is completely flat and vanishes for negative times 
(we estimate the scaled data has a relative error of less 
than 1%). The Bloch oscillations are always nearly as 
large as in the noninteracting case, and then they begin 
to decay. We are unable to determine whether they de- 
cay to a constant value as predicted by the Boltzmann 
equation, or whether there arc oscillations present in the 
steady state. In the quantum-mechanical system, there 
are two relevant time scales, the average time, and the 
relative time. The relative time governs the decay of 
the quasiparticle-like excitations, and this decay becomes 
rapid as the scattering increases. The average time gov- 
erns the Bloch oscillations, and it is not obvious from 
either the formalism or our results whether the steady- 
state current must be a constant, or whether it can os- 
cillate if the electric field is large enough (of course, with 
a period as short as these Bloch oscillations would have, 
they could not be measured in an experiment). 

In Fig- 01 we show results for the current with U = 0.5 
and U = 1. We take the time window to be larger here 
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FIG. 3: Nonequilibrium current for (a) U = 0.5 (At = 0.1) 
and (b) U = 1 (scaled from At = 0.1 and 0.0667) and longer 
times. Note how the Bloch oscillations are still present but 
become more erratic at large times. The current is multiplied 
by either 6 or 10 to enhance it at large times. 




Time 



FIG. 4: Nonequilibrium current for U = 2 and At = 0.0333. 
Note how the Bloch oscillations are now quite irregular and 
how there is substantial current before t = 0, because the 
data is not scaled. 
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to see if we can shed any further light on how the data 
evolves to the steady state but the time window is still too 
short. In Fig.0] we plot the current as a function of time 
for the small-gap Mott insulator with U — 2. It is much 
harder to achieve convergence for these results, and the 
scaling approach does not seem to work well, as the calcu- 
lated moments are less accurate for the scaled data, than 
for the data at the smallest At value (which is 0.0333 for 
the U = 2 data; relative errors here are probably at the 
10% level). Note how there is nonzero current at negative 
times, indicating that this data is not quite as accurate 
as the data with smaller U values. Note further, that 
the oscillatory behavior is quite irregular here, and it is 
not an exponentially decaying Bloch oscillation anymore. 
The Bloch oscillations rapidly change their character as 
the metal-insulator transition is crossed. 

Summary. We have developed the formalism for 
nonequilibrium dynamical mean-field theory. The basic 
method is similar to that of the equilibrium case, ex- 
cept we need to work in a real-time representation for 
all Green's functions, self-energies, and dynamical mean 
fields. The summation over momentum to yield the local 
Green's function now involves a two-dimensional integra- 
tion of a matrix-valued integrand. We presented numer- 
ical results for the Bloch oscillations in the presence of a 
large electric field, and showed how they decay as a func- 
tion of time when there is electron-electron scattering. 
In the transient response calculations presented here, we 
cannot determine whether oscillations are present in the 
steady state. We are currently working on a steady-state 
formalism should be able to shed light on that issue. 
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